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ABSTRACT 

A general purpose unstructured mesh solver for steady-state two-dimensional inviscid and 
viscous flows is described. The efficiency and accuracy of the method are enhanced by the 
simultaneous use of adaptive meshing and an unstructured multigrid technique. A method for 
generating highly stretched triangulations in regions of viscous flow is outlined, and a pro- 
cedure for implementing an algebraic turbulence model on unstructured meshes is described. 
Results are shown for external and internal inviscid flows and for turbulent viscous flow over a 
multi-element airfoil configuration. 
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1. INTRODUCTION 

Numerical methods for the computation of steady-state compressible flows have pro- 
gressed to the point where the major obstacle to achieving an efficient and accurate flow solu- 
tion for a given problem lies in one’s ability to generate an adequate mesh over the given 
geometry. For complex configurations, the most often employed approach consists of partition- 
ing the domain into a number of topologically simple regions and generating a structured mesh 
in each of these regions. The construction of these so-called block-structured meshes has pro- 
ven to be an expensive and time consuming process, requiring much human intervention. 

An alternate approach is afforded by the use of unstructured meshes using triangular ele- 
ments in two dimensions, and tetrahedral elements in three dimensions. Compared with struc- 
tured meshes, in which grid lines must propagate across the entire domain, unstructured meshes 
are purely local constructions, and thus provide a much greater degree of flexibility in meshing 
complex geometries. Furthermore, they provide a natural setting for the use of adaptive mesh- 
ing techniques, where new mesh points are added and the grid is restructured locally in regions 
where the flow gradients are large. 

However, even in the case of inviscid flows, the use of unstructured meshes has often 
been impeded by sometimes questionable accuracy and inefficient solvers. The efficiency of 
these solvers is hindered by the use of indirect addressing required by the random data sets, 
and by the use of relatively simple solution algorithms, due to the difficulties associated in con- 
structing implicit solvers on unstructured meshes, which require the inversion of large sparse 
matrices. While inefficiencies due to random data sets, which generally result in a reduction by 
a factor of three of the number of floating point operations per second (Mflop rate) achievable 
on present-day supercomputers, cannot be avoided, efficient solution algorithms requiring near 
optimal computational complexity may be devised. In the present work, this has been achieved 
through the use of a multigrid strategy [1]. Provided a consistent discretization of the govern- 
ing equations is employed, unstructured meshes can be shown to provide extremely accurate 
solutions through extensive use of adaptive meshing [2, 3, 4,5]. 

For viscous flows, unstructured meshes have seldom been employed, or in certain cases, 
have been used only in the inviscid regions of the flow, as part of a hybrid approach where 
structured meshes are employed in the regions of viscous flow [6,7]. The strong directionality 
of the gradients in the viscous flow regions, and the requirements they impose on the mesh 
generation procedure, as well as the frequent use of algebraic turbulence models, appears to 
have provided a strong deterrent against the use of fully unstructured meshes for viscous flows. 

In this paper, the development of an unstructured mesh solver for both inviscid and 
viscous steady-state flows about arbitrary two-dimensional geometries is described. This work 
represents an effort at constructing a general purpose, accurate, and efficient solution method. 
To this end, a general method for setting up the geometrical configuration (spline curves) and 
application of boundary conditions has been devised. An unstructured multigrid algorithm is 
used in conjunction with an adaptive meshing technique to ensure accurate and efficient solu- 
tions. For viscous flow calculations, a procedure for generating and adaptively refining highly 
stretched triangulations is presented. The implementation of an algebraic turbulence model for 
use on unstructured meshes is also described. 
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2. PROBLEM DEFINITION 

The first step involves the definition of the problem to be solved, which relates to the 
definition of the geometry, and the specification of boundary conditions and initial conditions. 

Geometry definition is the first step, which must necessarily be performed prior to the 
mesh generation and flow solution phases. In general, the geometrical configuration is 
described by an ordered series of points. These points are not employed as mesh points them- 
selves. A spline curve which fits these points is constructed, and mesh points are taken at 
predetermined locations along this spline curve. Thus, the geometrical configuration is in fact 
defined by a series of spline curves. This is especially important when adaptive meshing tech- 
niques are considered. New points which are introduced on the boundary must conform to the 
spline definition of this boundary, rather than simply being positioned midway between two 
neighboring grid points. Hence the geometry definition stage consists of identifying ordered 
groups of geometry points which are to be splined, selecting a type of spline for each group, 
and generating the spline curve for each group, which is stored as a series of spline coordinates 
and coefficients, for later use in the mesh generation and flow solution phases. 

The specification of boundary conditions forms part of the problem definition stage and 
as such, is best treated prior to the mesh generation and flow solution phases. The original 
edges which define the geometry (such edges link two geometry definition points prior to the 
spline fitting operation) are sorted into groups, each of which corresponds to a particular type 
of boundary condition. This information is then stored for subsequent use in the mesh genera- 
tion and flow solution phase. When the initial mesh is generated, boundary points are assigned 
the boundary condition corresponding to that of the geometry definition edge from which they 
were generated. When adaptively adding points to an existing mesh, new boundary points may 
also be assigned a boundary condition type in this manner. In the flow solution phase, applica- 
tion of the boundary conditions consists of looping through all boundary points and applying 
the appropriate boundary conditions at each point, as dictated by the boundary condition type 
associated with each point. To efficiently vectorize this step, boundary points are sorted into 
groups, each group representing a specific type of boundary condition, and the appropriate 
boundary condition is then executed on each group in a vector fashion. This predetermination 
and storage of boundaiy condition types allows for a completely general flow solution algo- 
rithm and facilitates the adaptive insertion of new boundary points without upsetting the logic 
of the solver. Finally the specification of initial conditions, such as Mach number, Reynolds 
number, and angle of incidence, may be effected without loss of generality as an input list in 
the flow solution phase. 

In this work, internal as well as external flow geometries have been considered. For 
external flow about multi-element airfoil configurations, each airfoil element is defined by a 
cubic spline. The circular far-field boundary is "splined" as a linear fit between a set of 
defining points. Tangential flow, or zero velocity boundary conditions are applied at the airfoil 
surfaces, depending on whether the Euler or Navier-Stokes equations are being solved. In the 
far field, a non-reflecting locally one-dimensional characteristic boundary condition is employed 
[ 8 ]. 

The internal flow about periodic cascade geometries requires the simultaneous use of two 
types of splines and three types of boundary conditions, thus providing a good illustration of 
the generality of this process. The initial definition of the geometry for this case, which is dep- 
icted in Figure 1, is composed of eight boundary regions. Regions 1 and 5 represent the 
inflow and outflow planes, regions 2,4 and 6,8 represent the periodic line. These regions are 
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thus all "splined" as straight line segments. Regions 3 and 7 represent the lower and upper sur- 
face of the turbine blade respectively, and are defined by cubic splines. In order to maintain 
continuity of slope and curvature at the leading and trailing edge points (i.e. the periodic points 
of the blade), in addition to the lower surface points, a periodic representation of the upper 
surface points is constructed and used in conjunction with these points to define the spline 
along the lower surface of the blade. A similar treatment is employed for the upper surface in 
boundary region 7. In regions 1 and 2, inlet and outlet boundary conditions are specified 
respectively. For the inlet flow, total pressure, total enthalpy, and the flow angle are specified, 
while the remaining condition is obtained by extrapolating the locally one-dimensional outgo- 
ing Riemann invariant normal to the boundary from the interior. For the outflow boundary, 
back pressure is specified, and total pressure, total enthalpy plus the outgoing locally one- 
dimensional Riemann invariant normal to the boundary are extrapolated from the interior 
[9,10]. In regions 3 and 7, flow tangency or zero velocity conditions are imposed, depending 
on whether inviscid or viscous flow solutions are sought. In regions 2, 4, 6, and 8, periodic 
boundary conditions are applied. A natural way of implementing periodic boundary conditions 
in the context of unstructured meshes is through the use of pointers. Each periodic boundary 
point is associated with its duplicate point on the corresponding periodic boundary through an 
integer pointer, which points to the address of the appropriate point. Thus, logically, the two 
corresponding periodic points refer to the same point, and are thus updated simultaneously and 
identically. However, duplicate physical copies of this point are required, each with a different 
physical coordinate, and each referring to a point on one of the two periodic cuts. 

3. MESH GENERATION 

3.1. Initial Mesh Generation 

The initial unstructured mesh is generated in three essentially independent stages. In the 
first stage, a distribution of mesh points filling the domain is generated. These points are then 
joined together to form a set of non-overlapping triangles using a Delaunay triangulation algo- 
rithm. Finally, a post-processing operation is employed to smooth out the mesh by slightly 
repositioning the points according to an elliptic smoothing operator [11]. 

While adaptive meshing techniques can be relied upon to increase the mesh resolution in 
regions of high flow gradients, a good initial mesh point distribution is essential to ensure the 
capture of all salient flow features on the initial mesh, and to reduce the number of adaptivity 
cycles required to attain a given accuracy level. This is particularly true in the case of high 
Reynolds number viscous flows, where very small normal spacings are required in the viscous 
regions. Since much effort has been expended in devising structured mesh generation strategies 
for specific types of geometries, these provide a natural starting point for the generation of a 
mesh point distribution. Each component of the geometry may be fitted locally with a struc- 
tured mesh, suitable to that particular type of component, and the union of all the points from 
these overlapping local meshes, which lie in the flow field, used as the basis for the triangula- 
tion. For multi-element airfoils, an O-mesh is fitted around each airfoil element using a hyper- 
bolic mesh generation algorithm [12]. For internal flow cascade geometries, the initial mesh 
point distribution is derived from a structured H-mesh [10]. 

Delaunay triangulation represents a unique way of joining a set of points in a plane 
together to form a set of non-overlapping triangles. Bowyer’s algorithm [13,14] is used to con- 
struct the triangulation. Assuming an initial triangulation exists (this may be constructed by 
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joining a small number of boundary points together), the mesh points are introduced and tri- 
angulated one at a time into the existing triangulation. Bowyer’s algorithm makes use of the 
circumcircle property of a Delaunay construction, which states that the circumcircles of the tri- 
angles may not contain vertices from other triangles. Thus, each time a new mesh point is 
introduced, the union of all triangles whose circumcircles contain this new point is identified. 
The existing mesh structure is removed in this region and new triangles are formed by joining 
the new point to all the vertices of the restructured region. When all mesh points have been 
introduced and triangulated, the initial unstructured mesh is obtained. 

3.2. Adaptive Mesh Generation 

This sequential insertion and local restructuring of new mesh points is ideally suited for 
an adaptive meshing strategy. Once a solution has been obtained on the initial mesh, new 
mesh points are created midway along mesh edges in regions of large flow gradients. Each 
new mesh point is then triangulated into the existing mesh using Bowyer’s algorithm. The 
search for the triangles whose circumcircles are intersected by the new point can be made 
extremely efficient by beginning with the two triangles on either side of the edge within which 
the new point was generated, and then searching through neighboring triangles. Hence, adap- 
tive mesh enrichment may be accomplished using only local searching and restructuring, thus 
avoiding the need for global mesh regeneration. When new boundary points are created, they 
must be positioned on the spline curve which defines that boundary. For concave boundaries, 
this results in mesh points which are not enclosed by any of the existing mesh triangles, as 
shown in Figure 2. To avoid failure of the intersected circumcircle search routine, new boun- 
dary points are initially positioned midway along the boundary edge within which they are 
generated. The circumcircle search and local retriangulation are then effected. The new mesh 
point is then displaced onto the spline curve, thus forming a sliver triangle which joins the new 
point with the two ends of the generating mesh edge. Such sliver triangles, which lead to a 
crossing of grid lines for concave boundaries (c.f. Figure 2), and for convex boundaries 
represent elements exterior to the computational domain, must be identified and subsequently 
removed. 


4. FLOW SOLUTION 

In conservative form, the full Navier-Stokes equations read 

dw | dfc dg c VyM, \ d/, 3s„l 

dt dx + dy Re_ [d* + J (*) 

where w is the solution vector and f c and g c are the cartesian components of the convective 
fluxes 

c f 

P P“ pv 

P u p M 2 + p pvu 

pv * c = p«v Sc= pv 2 + p (2) 

p£ J [ p uE+up J pvE+vp 

In the above equations, p represents the fluid density, u and v the x and y components of fluid 
velocity, E the total energy, and p is the pressure which can be calculated from the equation of 
state of a perfect gas 
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The viscous fluxes/, and g, are given by 
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where a represents the stress tensor, and q the heat flux vector, which are given by the consti- 
tutive equations for a Newtonian fluid 
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y is the ratio of specific heats of the fluid, M. the freestream Mach number, Re. the Reynolds 
number based on the airfoil chord, and Pr the Prandtl number. The coefficient of viscosity p 
varies with the temperature of the fluid, and is calculated as 

p = K r 0 72 ( fi ) 

where AT is a constant. Equation (1) represents a set of partial differential equations which 
must be discretized in space in order to obtain a set of coupled ordinary differential equations, 
which can then be integrated in time to obtain the steady-state solution. 

The spatial discretization procedure begins by storing flow variables at the vertices of the 
triangles. The stress tensor a and the heat flux vector q must be calculated at the centers of the 
triangles. This is achieved by computing the required first differences in the flow variables 
(from equations (5)) at the triangle centers. For a piecewise linear approximation of the flow 
variables in space, the first differences are constant over each triangle, and may be computed 

as 

1 „dw If 1 ^ Wfcrt + W k 

y ‘ A$ Wdy = 2 


-Cy* + i - ?*) 


(7) 


w, 




( 8 ) 

dy UMiy “ 4^ * ™ ~ 

where the summation over k refers to the three vertices of the triangle. The flux balance equa- 
tions are obtained by a Galerkin finite-element type formulation. The Navier-Stokes equanons 
are first rewritten in vector notation 

20L + v F = — — 

dt c Re. 

where the bold typeset denotes vector quantities. F c is a dyadic (second order tensor), the 


V.F, 


(9) 
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cartesian components of which are given by the/, and g, convective flux vectors deflned previ- 
ously. and a stmtlttr notation is employed for the viscous flux terms. Multiplying by a test 
function <(>, and integrating over physical space yields 


^j^wdxdy + J|<frV.F c dxdy = ^1=. fj$V.F v dxdy 
Integrating the flux integrals by parts, and neglecting boundary terms gives 




VvAf 

F c -V(J> dxdy - J|F V .Vd> dxdy 


( 10 ) 
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In order to evaluate the flux balance equations at a vertex P, <j> is taken as a piecewise linear 
function which has the value 1 at node P, and vanishes at all other vertices. Therefore, the 
integrals in the above equation are non-zero only over triangles which contain the vertex P 
thus defining the domain of influence of node P, as shown in Figure 3. To evaluate the above 
Integra s, we make use of the fact that <t> x and 4> y are constant over a triangle, and may be 
eva uated as per equations (7) and (8). The convective fluxes F c are taken as piecewise linear 
lunctions in space, and the viscous fluxes F v are piecewise constant over each triangle, since 

they are formed from first derivatives in the flow variables. Evaluating the flux integrals with 
these assumptions, one obtains h 




( 12 ) 

where the summations are over all triangles in the domain of influence, as shown in Figure 3 
da, "=P"*ents the directed (normal) edge length of the face of each triangle on the outer bourn 

tZ T' ■* * are ■“* « *• «> vertices *. eiflter end of Ws 

edge, and F. ts the v, scons flux in triangle e, e being a triangle in the domain of influence of * 

If the integral on the left hand side of equation (12) is evaluated in the same manner, the time 
totvanves become coupled in space. Since we are no. interested in the nme-accuracy ofte 

mmx' ^ s,ead )'; sl “' «* employ the concept of a lumped mass 

\ S 1S ec J ulvaIent to assuming w to be constant over the domain of influence while 
integrating the left hand side. Hence, we obtain h 

y/yM* 


d\V n n P? + X?B 

O L - V c c Af 

p dt J 2 ALab ~ 


Re. - 2 


E f (^ AL xii) 


(13) 

where die factor of 1/3 is introduced by the integration of * over the domain, and £1, represents 
die surface area of the domain of influence of P. For the convective fluxes, this procedure is 
quiva ent to the vertex finite-volume formulation described in [1,11 J. For a smoothly varying 
regular tnangulation, the above formulation is second-order accurate. S 

Additional artificial dissipation terms are required to ensure stability and to capture 
shocks without producing numerical oscillations. This is necessary for both inviscid and 
viscous flow computations, since in the later case, large regions of the flow field behave essen- 
tially mviscidly and the physical viscosity is not sufficient to guarantee numerical stability for 
e type of mesh spacings typically employed. Artificial dissipation terms are thus constructed 
as a blend of a Laplacian and a biharmonic operator in the conserved flow variables. The 
ap acian term represents a strong formally first-order accurate dissipation which is turned on 
only in the vicinity of a shock, and the biharmonic term represents a weaker second-order 
accurate dissipation which is employed in regions of smooth flow [11,15]. The spatially 
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discretized equations are integrated in time to obtain the steady-state solution using a five-stage 
time-stepping scheme, where the convective terms are evaluated at each stage within a time 
step, and the dissipative terms (both physical and artificial) are only evaluated at the first, third, 
and fifth stages. This particular scheme has been designed to maintain stability in regions 
where the flow is dominated by viscous effects, and to rapidly dampen out high-frequency 
error components, which is an essential feature for a scheme intended to drive a multigrid algo- 
rithm [16,17], Convergence is accelerated by making use of local time-stepping, implicit resi- 
dual averaging [16], and an unstructured multigrid algorithm [1]. 

The idea of a multigrid strategy is to accelerate the convergence to steady-state of a fine 
grid solution through corrections computed on coarser grids. An initial time step is performed 
on the fine grid, and the flow variables and residuals are then transferred to the coarse grid. A 
correction equation is constructed on the coarse grid by adding a forcing function to the origi- 
nal discretized equations. This forcing function is formed by taking the difference between the 
transferred residuals and the residuals of the transferred variables, thus ensuring that the evolu- 
tion of the coarse grid equations is driven by the fine grid residuals. Hence, when the fine grid 
residuals vanish, the coarse grid equations are identically satisfied, and generate zero correc- 
tions. After transferring values down from the fine grid, a time step is performed on the coarse 
grid, and the new values are transferred down to the next coarser grid. When the coarsest grid 
is reached, the computed corrections are successively interpolated back up to the finest grid, 
and the entire cycle is repeated. In the context of unstructured meshes, a sequence of coarse 
and fine meshes is best constructed by generating the individual meshes independently from 
one another (as opposed to subdividing a coarse mesh). Thus, in general, the coarse and fine 
meshes of a given sequence do not have any common mesh points or nested elements. Thus, 
the patterns for transferring the variables, residuals, and corrections back and forth between the 
various meshes of the sequence must be determined in a preprocessing operation, where an 
efficient tree-search algorithm is employed [1]. 

Such a multigrid algorithm may be combined with an adaptive meshing strategy in a 
natural manner. First, a sequence of globally generated meshes is constructed, and multigrid 
time-stepping is performed on this sequence until a satisfactorily converged solution is 
obtained. At this point, a new adaptively refined mesh is generated, and the transfer patterns 
for transferring variables from the previous mesh to the new mesh are determined. The flow 
variables are then transferred to this new mesh, providing a starting solution, and multigrid 
time-stepping is resumed on this new sequence which now contains an additional fine mesh. 
The process may be repeated, as shown in Figure 4, each time adding a new finer mesh to the 
sequence, until a converged solution of the desired accuracy is obtained. 

5. INVISCID FLOW RESULTS 

When the viscous terms in the Navier-Stokes equations are neglected, the Euler equations 
are obtained. Inviscid flow calculations can thus be computed using the method previously 
described, but neglecting the terms on the right-hand-side of equations (1) and/or (9), and 
replacing the no-slip wall boundary condition with a tangential slip velocity boundary condi- 
tion. Inviscid flow computations can be performed for a significantly lower cost than viscous 
flow computations, not only due to the reduced number of terms which need to be discretized, 
but also due to the elimination of the boundary layer and wake regions, where extremely high 
gradients generally occur, and which must be resolved in the viscous case. 
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5.1. External Flow Geometry 

The first test case consists of the computation of the inviscid compressible flow about a 
high-lift three-element airfoil configuration. The free-stream Mach number is 0.2, and the 
incidence is 8 degrees. At these conditions, the flow is entirely subcridcal. However, an 
extreme double suction peak occurs at the leading edge of the main airfoil, where the flow 
becomes nearly sonic. The capture of this extremely high localized gradient requires a very fine 
mesh resolution in this region. A sequence of seven meshes were used in the multigrid algo- 
rithm. The first three meshes were generated globally, and the further four meshes were gen- 
erated by adaptive refinement The criterion for adaptive mesh enrichment is based on the undi- 
vided difference of density [18]. The first difference of density computed along each mesh 
edge is examined. When this difference is larger than some fraction of the RMS average of all 
density differences across the mesh, a new point is added midway along the edge. A second 
pass is then performed which splits the remaining edges of each triangular element bordering 
on a previously split edge, thus ensuring an isotropic refinement. The finest mesh of this cal- 
culation is depicted in Figure 5. It contains 11949 nodes, of which 512 are on the airfoil sur- 
faces. Extreme refinement is seen to occur in the main airfoil leading edge region and in the 
gap regions. A globally refined mesh of this resolution would have required 10 to 20 times 
more points, and thus would be prohibitively expensive. The coarsest mesh of the sequence 
contains merely 110 points. The computed Mach contours in the flow field are depicted in Fig- 
ure 6, where the high gradients in the leading-edge region are evident. In Figure 7, a com- 
parison of the computed surface pressure distribution with that generated by a finite-element 
full potential solver [19], shows good agreement between the two methods, and illustrates the 
magnitude of the suction peak, where the pressure coefficient rapidly attains a value of -16.0. 
The convergence history for this case is depicted in Figure 8, where the fine grid residuals 
were reduced down to machine zero (in double precision) in just over 300 multigrid cycles. 
Each multigrid cycle required roughly 1.8 CPU seconds on a single processor of a Cray-2 
supercomputer, so that engineering calculations (50 - 100 mg cycles) could be obtained in 2 to 
3 CPU minutes. 

5.2. Internal Flow Calculations 

In a second test case, the inviscid flow through a turbine blade cascade geometry has 
been computed. The particular blade geometry has been the subject of an experimental and 
computational investigation at the occasion of a VKI lecture series [20]. A total of seven 
meshes were used in the multigrid algorithm, with the last three meshes generated adaptively, 
using the undivided density difference criterion. The coarsest mesh of the sequence contains 
only 51 points, while the finest mesh, depicted in Figure 9, contains 9362 points. Extensive 
mesh refinement can be seen to occur in the neighborhood of shocks, and in other regions of 
high gradients. The inlet flow incidence is 30 degrees, and the average inlet Mach number is 
0.27. The flow is turned 96 degrees by the blades, and the average exit isentropic Mach 
number is 1.3. At these conditions, the flow becomes supersonic as it passes through the cas- 
cade, and a complex oblique shock wave pattern is formed. These are evident from the com- 
puted Mach contours depicted in Figure 10. All shocks are well resolved, including some of 
the weaker reflected shocks, which non-adapted mesh computations often have difficulty 
resolving [10]. Details of the flow in the rounded trailing edge region of the blade are shown 
in Figure 10, where oblique shock waves are formed on the upper and lower surfaces of the 
blade, and where the flow separates (inviscidly), forming a small recirculation region. The 
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surface isentropic Mach number distribution for this case is compared with experimental data 
provided from [20], in Figure 1 1 . The upper surface experimental values correspond to an exit 
isentropic Mach number of 1.31, while the lower surface experimental values are taken from a 
case where the measured exit isentropic Mach number was 1.21. These values are included 
*ince lower surface data was not available for an exit Mach number of 1.31, and since the 
lower surface values were seen to be relatively insensitive to the exit Mach number in this 
range. Keeping in mind the inviscid nature of the computational results, good agreement is 
observed over most of the surface of the blade. The irregularities in the numerical solution are 
due to a poor surface definition of the blade, the effect of which is enhanced by the adaptive 
meshing procedure, which tends to refine the mesh in the vicinity of non-smooth geometries. 
Once the first four globally generated meshes were constructed, the entire flow solution - adap- 
tive mesh enrichment cycle was performed three times, executing 25 multigrid cycles at each 
stage. This entire operation required 40 CPU seconds on a single processor of a Cray-YMP 
supercomputer. The residuals on the finest mesh were reduced by two and a half orders of 
magnitude, which should be adequate for engineering calculations. The efficiency of this solu- 
tion illustrates the possibility of constructing a truly interactive adaptive mesh Euler solver in 
two-dimensions, provided the present algorithm may be efficiently implemented in parallel on 
all eight processors of the Cray-YMP. 

6. CONSIDERATIONS FOR HIGH REYNOLDS NUMBER VISCOUS FLOWS 

While the methodology previously described applies in principle to viscous flows as well, 
the efficient solution of high Reynolds number flows requires the generation of highly stretched 
meshes as well as the implementation of a turbulence model for use on unstructured meshes. 

6.1. Mesh Generation 

The generation of highly stretched unstructured meshes requires a suitable mesh point 
distribution, with closely packed points in the normal direction, and sparsely distributed points 
in the streamwise direction, as well as a method for producing an appropriate triangulation of 
such a point distribution. The generation of a suitable point distribution may be effected as pre- 
viously described, using local hyperbolically generated structured meshes. However, a 
Delaunay triangulation of a given set of points tends to produce the most equiangular triangles 
possible, and therefore in general, is not well suited for the generation of highly stretched mesh 
elements. Thus, an alternate triangulation procedure must be employed. The approach taken 
consists of defining a stretching vector (stretching magnitude and direction) at each node of the 
initial point distribution throughout the flow field. Assuming an initial triangulation has been 
obtained, when a new mesh point is to be inserted, the associated stretching vector is employed 
to construct a locally mapped space such that, within this mapped space, the local point distri- 
bution appears isotropic. A Delaunay triangulation is then performed to triangulate the new 
point into the mesh in this mapped space, and the resulting triangulation is mapped back into 
physical space, thus resulting in the desired stretched triangulation [21], Hence, a fully 
unstructured mesh with highly stretched elements in the boundary layer and wake regions, 
nearly equilateral triangles in the inviscid regions of flow, and a smooth variation of elements 
throughout the transition regions is obtained. The use of fully unstructured meshes for viscous 
flow calculations has been pursued, as opposed to the hybrid structured-unstructured meshes 
often advocated in the literature [6,7], due to the increased generality they afford in dealing 
with geometries with close tolerances between neighboring bodies, where confluent boundary 



- 10 - 


layers may occur, and due to the ease with which adaptive meshing may be incorporated 
throughout the viscous and inviscid regions of flow. 

6 . 2 . Turbulence Modeling for Unstructured Meshes 

Algebraic turbulence models typically require information concerning the distance of each 
mesh point from the nearest wall. Turbulence length scales, which are related to the local 
boundary layer or wake thickness, are determined by scanning the appropriate flow values 
along specified streamwise stations. For example, the Baldwin-Lomax turbulence model [22] 
uses the location of the maximum of the moment of vorticity along streamwise stations normal 
to the boundary layer to estimate the local turbulence length scales. In the context of unstruc- 
tured meshes, mesh points and thus flow variables do not naturally occur at regular streamwise 
locations. Thus, lines normal to the walls and viscous layers must be created and flow variables 
interpolated onto these lines, in order that the turbulence length scales may be determined. This 
type of approach has previously been implemented for supersonic ramp geometries by Rostand 
[23]. However, in the present work, more complex geometries must be accommodated. Recal- 
ling that, in the mesh generation procedure, the initial mesh point distribution was obtained by 
generating a series of local structured meshes about each geometry component, a natural 
manner of creating streamwise turbulence modeling stations is to make use of these local struc- 
tured meshes as background meshes for the turbulence modeling routine. Thus, at each time 
step in the solution procedure, the current flow variables from the global unstructured mesh are 
interpolated onto each of the local structured meshes. The Baldwin-Lomax turbulence model is 
executed on these local structured meshes, and the resulting eddy viscosity distribution is inter- 
polated back onto the unstructured mesh. In regions where two or more local structured meshes 
overlap, the multiple eddy viscosity values (one from each local mesh) which are interpolated 
back to the unstructured mesh are each weighted by their relative distance from the respective 
wall. Structured mesh lines emanating from one geometry component are terminated if they 
intersect a neighboring component, as shown in Figure 12, such that in any region of the flow 
field, the eddy viscosity is only related to the viscous layers and associated walls which are 
directly visible from that point. The same interpolation routines employed for the multigrid 
algorithm are used to pass variables back and forth between the global unstructured mesh and 
the local structured turbulence meshes. The determination of these transfer patterns is done in 
an efficient preprocessing operation, and the transfer addresses and coefficients are stored for 
subsequent use in the turbulence modeling routine. The whole process is very efficient, and in 
general, the entire turbulence modeling routine, including the interpolation procedures requires 
only 15% of the total time within a multigrid cycle. Memory requirements are however 
increased by about 50%, since extra variables and transfer coefficients must be stored for the 
local structured meshes. 

7. TURBULENT FLOW RESULTS 

The above modifications have been incorporated into the present scheme in order to com- 
pute the turbulent flow in the transonic regime past a two-element airfoil. The configuration 
consists of a main airfoil with a leading edge slat, which has been the subject of extensive 
wind tunnel tests [24], as part of a program aimed at improving the maneuvering capabilities of 
fighter aircraft in the transonic regime. A multigrid sequence of five meshes was employed to 
compute the flow about this configuration. The finest mesh is depicted in Figure 13. It con- 
tains 22,509 points, of which 256 lie on the surface of the main airfoil, and 128 on the surface 
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of the slat. The average width of the elements on the airfoil surfaces is 0.00001 chords, result- 
ing in cell aspect ratios of the order of 1000:1 in these regions. The turbulence background 
meshes, consisting of one local structured mesh for each airfoil component, are depicted in 
Figure 12. The computed Mach contours are shown in Figure 14. For this case, the ffeestream 
Mach number is 0.5, the Reynolds number based on the chord is 4.5 million, and the incidence 
is 2.8 degrees. At these conditions, the flow is supercritical, and a shock is formed on the 
upper surface of the slat, as can be seen from the figure. The computed shock is somewhat 
diffuse, and an increased mesh resolution in this region would be required to obtain a crisper 
shock definition. However, the sudden thickening of the boundary layer as it interacts with the 
shock is evident from the computed Mach contours. A small recirculation region is also 
observed on the lower surface of the slat. A comparison of the computed surface pressure dis- 
tribution with the experimental wind tunnel data is given in Figure 15. Computed and experi- 
mental values are seen to agree favorably in all regions, demonstrating a good prediction of the 
suction peaks, and location of the slat upper surface shock. This solution required roughly eight 
cpu minutes on a single processor of a Cray-2 supercomputer, which corresponds to 75 mul- 
tigrid cycles on the finest grid, during which the residuals were reduced by approximately three 
orders of magnitude. To the author’s knowledge, this represents the first compressible turbulent 
flow calculation for multi-element airfoil geometries using unstructured meshes. 

8. CONCLUSION 

A method for solving viscous and inviscid flows about arbitrary two-dimensional 
configurations has been presented. An attempt has been made to keep the method as general as 
possible. The simultaneous use of multigrid and adaptive meshing results in a rapidly conver- 
gent and accurate solution. For a given number of unknowns (mesh points), unstructured mesh 
solutions can be obtained with roughly the same number of operations as is required by the 
most efficient current structured mesh solvers. However, the speed of execution of unstruc- 
tured mesh codes on present-day vector computers is roughly three times slower than that 
observed with structured mesh codes, due to the indirect addressing and scatter-gather opera- 
tions required by the use of random data-sets. However, this factor can easily be outweighed 
through the use of a more efficient placement of grid points, using adaptive meshing tech- 
niques. For inviscid computations, the present algorithm appears robust and efficient enough 
that it may be implemented in an interactive mode on the latest generation of supercomputers. 
An inviscid solver based on these techniques in three dimensions, which is in the planning 
stages, should also provide a competitive solution technique for large problems. For turbulent 
viscous flow calculations, substantial modifications to the mesh generation phase were required. 
The implementation of an algebraic turbulence model has demonstrated a good prediction capa- 
bility for flow over streamlined bodies. In future work, the implementation of a more general 
turbulence model, such as a field equation model, will be considered for flow over arbitrary 
geometries with massive separation. 
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Figure 3 

Domain of Influence of Finite-Element Basis Function and Equivalent 
Finite-Volume Control Volume 
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Figure 8 

Convergence Rate on Finest Mesh for the Three-Element Airfoil Case 
as Measured by the Average of the Density Residuals throughout the Flowfield 
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Figure 10 

Computed Mach Contours for Flow Through a Periodic Turbine Blade 

Cascade Geometry 
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Figure 11 

Comparison of Computed Surface Isentropic Mach Number Distribution 
with Experimental Values for Flow Through a 
Periodic Turbine Blade Cascade Geometry 
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Flgure 12 

Illustration of Local Structured Background Turbulence Modeling Meshes Employed 
for Computing Turbulent Flow over a Two-Element Airfoil Configuration 
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Figure 14 

Computed Mach Contours of Supercritical Turbulent Viscous Flow over a Two-Element Airfoil Configurat.on 
Mach Number = 0.5, Reynolds Number = 4.5 million. Incidence = 2.8 degrees 
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Figure 15 

Compjmson of Computed Surface Pressure Distribution with Experimental 
Wind-Tunnel Data for Flow Over Two-Element Airfoil Configuration 
Mach Number = 0.5, Reynolds Number = 4.5 million, Incidence = 2.8 degrees 
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